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The planar Poiseuille flow induced by a constant external field (e.g., gravity) has been the subject 
of recent interest in the case of molecular gases. One of the predictions from kinetic theory (confirmed 
by computer simulations) has been that the temperature profile exhibits a bimodal shape with a 
local minimum in the middle of the slab surrounded by two symmetric maxima, in contrast to the 
unimodal shape expected from the Navier-Stokes (NS) equations. However, from a practical point 
of view, the interest of this non-Newtonian behavior in molecular gases is rather academic since 
it requires values of gravity extremely higher than the terrestrial one. On the other hand, gravity 
plays a relevant role in the case of granular gases due to the mesoscopic nature of the grains. In 
this paper we consider a dilute gas of inelastic hard spheres enclosed in a slab under the action of 
gravity along the longitudinal direction. In addition, the gas is subject to a white-noise stochastic 
force that mimics the effect of external vibrations customarily used in experiments to compensate 
^ . for the collisional cooling. The system is described by means of a kinetic model of the inelastic 

Boltzmann equation and its steady-state solution is derived through second order in gravity. This 
solution differs from the NS description in that the hydrostatic pressure is not uniform, normal 
stress differences are present, a component of the heat flux normal to the thermal gradient exists, 
and the temperature profile includes a positive quadratic term. As in the elastic case, this new term 
is responsible for the bimodal shape of the temperature profile. The results show that, except for 
high inelasticities, the effect of inelasticity on the profiles is to slightly decrease the quantitative 
deviations from the NS results. 
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As is well known, the Poiseuille flow is a typical example of fluid dynamics described in many textbooks. 1 In 
£S) | its classical formulation, the Poiseuille problem consists of finding the flow velocity and temperature profiles of a 
^ ■ Newtonian fluid enclosed in a slab or in a pipe and subject to a longitudinal pressure gradient. Essentially the same 
effect is generated when the longitudinal pressure difference is replaced by a uniform gravitational force mg directed 
On : longitudinally. This latter mechanism for driving the Poiseuille flow does not produce longitudinal gradients and so 
has proven to be more convenient than the former in computer simulations as well as from the theoretical point of 
view, especially to assess the reliability of the continuum d e scriptioni2i2ii£i£iL2i2ii2iiiiI2ii&I2ii£iI& 

Kinetic theory analyses of the gravity-driven Poiseuille flow based on an expansion in powers of the gravity strength 
g^gilSjl^ on Grad's moment methodf2iii or on an expansion in powers of the Knudsen number, ^ show interesting 
"j***; i non-Newtonian effects. In particular, to second order in g the temperature profile includes a positive quadratic term, 
in addition to the negative quartic term predicted by the Navier-Stokes (NS) description. As a consequence of this new 
term, the temperature profile exhibits a bimodal shape with a local minimum at the middle of the channel surrounded 
by two symmetric maxima at a distance of a few mean free paths. In contrast, the NS hydrodynamic equations 
predict a temperature profile with a (flat) maximum at the middle. The Fourier law is dramatically violated since in 
the slab enclosed by the two maxima the transverse component of the heat flux is parallel (rather than anti-parallel) 
, to the thermal gradient. This correction to the NS temperature profile is not captured by the next hydrodynamic 
description, i.e., by the Burnett equations^S The kinetic theory prediction of a bimodal temperature profile has 
been confirmed by numerical Monte Carlo simulations of the Boltzmann equationiiS*i^ and by molecular dynamics 
simulations On the other hand, when the Poiseuille flow is driven by a longitudinal pressure gradient instead of an 
external force, the NS description is in good agreement with Monte Carlo simulations of the Boltzmann equation^ 
Notwithstanding its theoretical and academic interest, the Poiseuille flow induced by gravity is of little practical 
interest for conventional gases under terrestrial conditions. At a microscopic level, the relevant dimensionless param- 
eter measuring the strength of gravity is g\/v% h , where A is the mean free path and Vth is a typical molecular speed 
(or thermal velocity). The parameter g\/v^ h measures the effect of gravity on a molecule between two successive col- 
lisions. For instance, in the case of argon at room pressure and temperature, one has A ~ 700 A and Vth ~ 400 m/s^i 
so that g\/vl h ~5x 10~ 12 . 

The negligible effect of gravity on molecular gases is a consequence of their small mean free paths and large 
thermal velocities over mesoscopic or macroscopic scales. However, this is not necessarily so when dealing with a 
"granular" g aSf igii2i22i2ii22i2ii i.e., a collection of a large number of discrete solid particles (or grains) in a fluidized 
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state such that each particle moves freely and independently of the rest, except for the occurrence of inelastic binary 
collisions. Depending on the material properties of the grains, the solid fraction, and the state of agitation, the 
parameter g\/v% h can take values within a wide spectrum. Let us take three representative examples. In Ref. |2J, the 
statistical properties of stainless-steel spheres of diameter a — 3.175 mm rolling on an inclined surface and driven by 
an oscillating wall were experimentally studied. Typical values of the mean free path and the thermal velocity were 
A ~ 1 cm and vth ~ 1 cm/s, which leads to g\/v% h ~ 10 3 . Experiments on glass beads of diameter a = 4 mm driven 
by a vertically oscillating boundary were reported in Ref. 25. In those experiments, A ~ a and vth ~ 20 cm/s, so 
that g\/v^ h ~ 1. As a final example, we consider the experiments carried out in a flying rocket on bronze spheres of 
diameter a — 0.3-0.4 mm excited by vibrations^ From the experimental data corresponding to the most dilute cell 
one can estimate A = 1 mm and u t j, = 5 m/s; under terrestrial conditions (g — 9.8 m/s 2 ), this implies g\/v^ h ~ 10 . 

In this paper we address the granular Poiseuille flow generated by gravity under the assumption that g\/v^ h is (i) 
large enough as to produce noticeable gradients of density, flow velocity, and granular temperature, but (ii) small 
enough as to allow for a perturbative treatment; roughly speaking, this corresponds to 10~ 3 < gX/v^ h < 10 _1 . Since 
kinetic energy is continuously being dissipated by inelastic collisions, we assume that the gas is externally excited by 
a "heating" mechanism. This guarantees that the gas is in a (uniform) steady state even in the absence of gravity. 
As the simplest way of mimicking energy input through boundary vibrations, we consider the widely used stochastic 
force with white noise properties. This means that every particle receives uncorrelated random kicks. Besides, the 
relative magnitude of the kicks scales with the square root of the local collision rate. 

Our main goal is to derive the profiles of the hydrodynamic variables and their fluxes in the bulk region, and assess 
to what extent they are influenced by the degree of inelasticity. In principle, an adequate framework to undertake 
this task is provided by the Boltzmann equation for inelastic hard spheres. However, its mathematical intricacy 
prevents one from deriving practical results, even in the elastic case, unless Grad's method with a high number of 
momentsii or the direct simulation Monte Carlo methodic are employed. In order to get explicit expressions with 
a moderate calculation effort, we replace the Boltzmann inelastic collision operator by a much more tractable kinetic 
model recently proposed^! as an extension to granular gases of the celebrated Bhatnagar-Gross-Krook (BGK) model 
for conventional gases^ The resulting kinetic equation is solved through second order in g and the associated profiles 
of the hydrodynamic fields and their fluxes are derived. The results show that the same type of non-Newtonian 
properties that appear in the elastic case are present for granular gases as well. On the other hand, for small and 
moderate inelasticities, we observe that those effects tend to be slightly inhibited as the inelasticity increases. 

The organization of the paper is as follows. Section [H] is devoted to the description of the flow under study and 
its solution in an NS hydrodynamic description. The kinetic theory description is presented in Sec. IIIII where a 
perturbation expansion in powers of gravity is carried out. The results are summarized and discussed in Sec. IIVI 
Finally, the main conclusions of the paper are briefly presented in Sec. [V] 



II. STATEMENT OF THE PROBLEM 



A. Inelastic hard spheres 



Let us consider a granular gas composed of smooth inelastic hard spheres of diameter a, mass m, and coefficient of 
normal restitution a. In the dilute regime, the one-particle velocity distribution function /(r, v; t) obeys the (inelastic) 
Boltzmann equation22i22i 

4+v-V + gi+f)/= J[/,/], (2.1) 

where g is the acceleration due to an external force, T is the operator representing the action of a given heating 
mechanism to compensate for the collisional energy loss, and J[f, f] is the Boltzmann collision operator. Its expression 
is 

J[f, f] = a 2 J dvi J da 9(v 01 • <t)(v i ■ a) [a- 2 /(v")/(v'/) - /(v)/(vi)] , (2.2) 

where the explicit dependence of / on r and t has been omitted. In Eq. (|2.2(l . O is the Heaviside step function, a is 
a unit vector directed along the centers of the two colliding spheres at contact, vqi = v — vi is the relative velocity, 
and the pre-collisional or restituting velocities v" and v" are given by 



v" = v - 1 J~ a (vqi • a)a, w'l = vi + 1 J" 01 (vox ■ a)a. (2.3) 

2a 2a 
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The first few moments of the distribution function define the number density n, the flow velocity u, and the granular 
temperature T: 



n(r,t) \ - / 1 
n(r, t)u(r, t) = / dv i v 
n(r,t)T(r,t) J J \f\ 



v /(r,v;t), (2.4) 

m ^2 

where V = v — u is the velocity relative to the local flow. The macroscopic balance equations for the local densities 
of mass, momentum, and energy follow directly from Eq. 12. If) by taking velocity moments: 

D t n + nV- u = 0, (2.5) 



D t u+ V • P = g, 

mn 



D t T + — (V • q + P : Vu) = -(C - l)T. 



In these equations, D t = dt + u ■ V is the material time derivative, 

P(r,t)=m y dvVV/(r,v;t) 



is the pressure or stress tensor, 



is the heat flux, 



q(r,f) = f / rfvy 2 V/(r,v;i) 



C(r,t) 



?7i 



3n(r,t)T(r,t) 

is the cooling rate associated with the inelasticity of collisions, and 

m 



dW 2 J[f, /] 



7(r,*) = - 



3n(r,t)T(r 



J rfvy 2 f/(r,v;i) 



(2.6) 
(2.7) 

(2.8) 
(2.9) 
(2.10) 
(2.11) 



is the heating rate associated with the external driving J- . Upon writing Eqs. I|2.5fl and l|2.6|l it has been assumed 
that T preserves the local number and momentum densities, i.e., 



J dv^/(r,v;t) = J dvvf/(r,v;t) = 0. 



(2.12) 



Equation l|2.10|l shows that the cooling rate is a complicated nonlinear functional of /. By dimensional analysis, 
C oc nT 1 / 2 , but the proportionality constant is an unknown function of a. A reasonable estimate of ( can be obtained 
by replacing in Eq. H2.10|) the actual velocity distribution function / by its local Maxwellian approximation 



fi(r,v;t) = n(r,t) 



771 


3/2 


exp 


2nT(r,t)_ 



to (v — u(r, t)Y 



The result is^i 



where 



16 2 /7rT\ 
V = 77fT I I 



1/2 



5 \ m J 

is an effective collision frequency, independent of the coefficient of restitution a. 



(2.13) 



(2.14) 



(2.15) 
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B. Gravity-driven Poiseuille flow 



Now we assume that the granular gas is enclosed between two infinite parallel plates normal to the y-axis. A 
constant external force per unit mass (e.g., gravity) g = — gz is applied along a direction z parallel to the plates. The 
geometry of the problem is sketched in Fig. ^ 

As done in laboratory experiments (and in computer simulations), we will assume that energy is externally injected 
into the system to compensate for the collisional cooling, so that a steady state is achieved even if the gravity field 
is formally switched off. In real experiments; 2 ^**^ this is usually achieved by means of boundary vibrations of 
small amplitude A ~ a and high frequency ui/2tt ~ 10-100 Hz, so that the maximum accelaration T — Alu 2 is 
usually several times larger than the acceleration due to gravity on Earth. However, this type of realistic heating 
through the boundaries is difficult to deal with at a theoretical level due to unavoidable boundary effects. These 
difficulties are overcome by assuming a bulk heating mechanism acting on all the particles simultaneously. The most 
commonly used type of bulk driving for inelastic particles consists of a stochastic force in the form of Gaussian white 
noise22i2&i2!iSi2ii2L22iSi. More precisely, each particle i is subject to the action of a stochastic force Fj(£) that has the 
properties 

(F,(t)) = 0, (F i (t)F J (t')> = lm 2 £ 2 %<5(i - 1'), (2.16) 

where I is the 3x3 unit matrix and £ 2 represents the strenght of the correlation. According to this white noise driving, 
during a small time step St each particle i receives an independent "kick" such that its velocity is incremented by a 
random value 5vj with the statistical properties^ 

(Svt) = 0, {SvtSvj) = \fStSij. (2.17) 

Therefore, \Sv\ ~ £v6~i. The associated operator T appearing in the Boltzmann equation l|2.1[) is^i 



2\dv) 



(2.18) 
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Thus £ 2 /2 plays the role of a diffusion coefficient in velocity space. The operator (|2.f 8[> verifies the properties i|2.I 2fl . 
while insertion into Eq. (|2.f f|) shows that the heating rate is 

7=^- (2.19) 

It still remains to define the spatial dependence of 7. By simplicity we assume that the white noise driving 
compensates locally for the collisional energy loss. This means that 7 = £ or, equivalently, £ = W (T/m at any point. 
This choice can be justified by the following argument. Since, as seen above, |<5v| ~ £V~St, the choice 7 = C implies 
that 

^ ~ Jis5t(l-a 2 ), (2.20) 

Vth 



where use has been made of Eq. H2.I4J) and of Vth ~ y/T/m. Equation (|2.20|l means that the relative random increment 
of velocity at a given point scales as the square root of the average collision number at that point. When heating the 
gas through the boundaries, the energy input is propagated to the whole system by means of collisions. Since the 
white noise driving intends to mimic that effect, it is quite natural that the relative magnitude of the kicks is larger 
in those regions where the collisions are more frequent. 

By considering the above white noise excitation mechanism, a steady state can be expected in which the physical 
quantities depend on the coordinate y only and the flow velocity is parallel to the z axis, u = u z (y)z. In that case, 
the Boltzmann equation l|2.I|) becomes 

CT d 2 d d\ 

L^-^+^) /=j ^ (2 - 2i) 

Similarly, the balance equations for momentum and energy, Eqs. (|2.6|) and l|2.7|l . reduce to 

BP 

vv - 0, (2.22) 



pg, (2.23) 



dy 

9P V z . 
dy 

where p = ran is the mass density. Note that the inelasticity does not appear explicitly in the balance equations 

C. Navier— Stokes description 

In the Newtonian description the fluxes are related to the hydrodynamic gradients by the Navier-Stokes (NS) 
constitutive equations. 3140-41 In the geometry of the Poiscuille problem they read 

Pxx — Pyy = Pzz = Pj (2.25) 



dT dn 

q v = - K dy-- f %> (2 - 2?) 



qz -0, 



(2.28) 
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where p = nT = ^Tr P is the hydrostatic pressure, rj is the the shear viscosity, k is the thermal conductivity, and /x 
is a transport coefficient with no analog for elastic fluids. These transport coefficients can be explicitly derived from 
the Boltzmann equation (|2.1(l by application of the Chapman-Enskog method in the first Sonine approximation. In 
the case of the white noise heating lf2~T%j) with £ = ^C,T/m their expressions art 



.41 



'V 

where 



»? = -, «=^-(l + 2fc), V=^—k, (2.29) 
v„ lmv K Zmv K 



v , , / 49 - 33a 19 - 3a , . 
"« = o (1 + «) — T^" + ETo ■ (2-31) 



3 V 'V 16 512 

In the above equations, v is the effective collision frequency defined by Eq. I|2.15|l and k is the kurtosis of the 
homogeneous heated state. Its expression is well approximated by2£ 



16(1 - a)(l - 2a 2 ) 
241 - 177a + 30a 2 (1 - a)' 



(2.32) 



The kurtosis k is rather small for all a. In particular, \k\ < 0.013 for 0.6 < a < 1. Therefore, one can neglect k in 
(|2~2^|) - (|23T|) to get 

P 4 _5p 48 

77 ~ !/ (1 + a) (3 - a) ' K ~ 2m^ (1 + a) (49 - 33a) ' M " 1 j 

In the interval 0.6 < a < 1, the expressions l|2.33|l for 7/ and k deviate from those of (|2.29|) less than 0.04% and 3%, 
respectively. Besides, the ratio nfi/Tn is smaller than 0.013, so that \i can be neglected. Note that the negligible 
role played by /j does not hold in the homogeneous cooling state^i*^ It is worth pointing out that, while the 
shear viscosity monotonically increases with inelasticity, the thermal conductivity starts decreasing with increasing 
inelasticity, reaches a minimum value around a = 0.4, and then slightly increases for a > 0.4. This non-monotonic 
behavior of k in the heated state contrasts with the one found in the free cooling case^i^i*^ 
Combining Eqs. fTZfy - fTTfy . we get 

|=0, (2.34) 

3T>7r = «, (2 ' 35) 

dy dy 

°^ = -J^)\ (2 . 36 ) 



dy dy \ dy 

where k' = n — rifi/T ~ k. Equation l|2.35|l gives a parabolic- like velocity profile, that is characteristic of the Poiscuillc 
flow. The temperature profile has, according to Eq. H2.36fl . a quartic-like shape. Strictly speaking, these NS profiles 
are more complicated than just polynomials due to the temperature dependence of the transport coefficients. Since 
the hydrodynamic profiles must be symmetric with respect to the middle plane y — 0, their odd derivatives must 
vanish at y = 0. Thus, from Eqs. H2.35JI and (|2.36(l we have 



d 2 u z 



dy 2 



= Po9 d^T 
y =a Vo ' dy 2 



= 0, 

y=o 



d 4 T 



y=0 ^0^0 



2^L, (2.37) 



where the subscript denotes quantities evaluated at y — 0. According to Eq. (|2.37|) . the NS equations predict that 
the temperature has a maximum at the middle layer y = 0. As we will see in Sec. IIIII the kinetic theory description 
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shows that the temperature actually exhibits a local minimum at y — 0, since d 2 T/dy 2 \ y0 is a positive quantity (of 
order g 2 ). 

The closed set of nonlinear equations (|2.34|I - H2.36|I cannot be solved analytically for arbitrary g due to the spatial 
dependence of the transport coefficients. On the other hand, if the acceleration of gravity is sufficiently small at 
the microscopic scale, we can expand in powers of g and keep the first few terms only. To second order, the NS 
hydrodynamic profiles near the layer y — are 

u z (y) = u + ^y 2 + O(g 3 ), (2.38) 

T(y)=T --^ T y i + O(g i ). (2.39) 
1 2?7oKo 

The space variable y can be eliminated between Eqs. (|2.38|l and l|2.39|) to obtain the following nonequilibrium "equation 
of state" : 

T = T Q -^ T (u -u z ) 2 + O(g 4 ). (2.40) 

The NS profiles for the fluxes are 

P yz (y) = -Pogy + O( g 3 ), (2.41) 



q y (y)^^y 3 + 0( 9 i ). (2.42) 

III. KINETIC THEORY DESCRIPTION 
A. A kinetic model 

In this Section we will see that most of the NS predictions discussed in the preceding Subsection do not hold true, 
even to first order in g, when the problem is attacked from a more detailed kinetic point of view. In principle, the 
task consists of solving the Boltzmann equation 1)2. 21|) through order g 2 in a region near the central layer y = 0. 

Given the mathematical complexity of the Boltzmann collision operator (|2.2|l . especially in the case of inelastic 
collisions, we simplify the analysis by replacing J[/, /] by a BGK-like kinetic modeliSL^ 

J[f, f] - -(3(a)v(f - ft) + | A • [(v - u) /] , (3.1) 

where v is the collision frequency (|2.15|) , fi is the local Maxwellian distribution (|2.13|l , and Q is the associated cooling 
rate (|2.14|l . In addition, f3(a) is a dimensionless function of the coefficient of restitution that can be freely chosen to 
optimize agreement with the Boltzmann description. Equation (|3.1() differs from the original formulation of the model 
kinetic equation^! in that the exact (local) homogeneous cooling state of the Boltzmann equation is replaced by ft 
and the exact cooling rate 1)2. 10|) is approximated by Q. Confirmation of the quantitative agreement between the 
kinetic model and the Boltzmann equation has been found for the simple shear flo w 44 ' 45 and the nonlinear Couette 
flowi£ 

The first term on the right-hand side of l|3.1|l describes collisional relaxation towards the local Maxwellian with a 
collision rate (3v, while the second term describes the dominant collisional cooling effects. The necessity for this term 
to accurately represent the spectrum of the Boltzmann collision operator is discussed in Ref. OOJ. However, it can be 
viewed more simply as an effective "drag" force that produces the same energy loss rate as that produced by the 
inelastic collisions. The NS transport coefficients derived from the model (|3.1|l in the case of white noise heating are— 

5P M = 0. (3.2) 
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A simple choice for the parameter /3 is (3 — |(1 + a)m* On the other hand, comparison with the (approximate) 
Boltzmann results l|2.33[l shows that the shear viscosity is reproduced if j3 takes the value 

2 + Q 



0= (1 + a)- 



(> 



while the thermal conductivity is reproduced if 



P = (1 + a) 



19 - 3a 



(3.3) 



(3.4) 



The discrepancy between Eqs. (|3.3|l and (|3.4|) persists in the elastic limit (a = 1) and is a well-known limitation of 
the BGK model. As will be seen in Sec. IIVI one can partially circumvent this problem by expressing the final results 
in terms of 77 and k. 

Inserting the model (|3.1|) into Eq. I|2.21[l . we get the kinetic equation 



d d\ 



m ov , 



where, for consistency, we have made the approximation £ 
from the local equilibrium distribution, let us write 



V 1 /'. (3.5) 

Q in Eq. 12.21[1 . In order to focus on the deviations 



/ = .Mi + $). 



(3.6) 



Then, Eq. I|3.5[) becomes 



(1 + $) 



du 2 



Vyd y log f t - [ g + V y — ) — log ft 



0_ 

dVz 



■ + v« 



di 



d 



dy J dV z 



4> 



Vydy* 



9V 



(3.7) 



where the operator d y = d/dy + (du z /dy)d/dV z derives with respect to y at constant V (i.e., not at constant v). 
Moreover, in Eq. (|3.7|) we have introduced the modified collision frequency v' = (3h> + Q. As Eq. (|3.2|) shows, v' is the 
effective collision frequency associated with the shear viscosity of the heated granular gas in the kinetic model. 

Since we are interested in the solution of Eq. 1)3. 7JI in the bulk, it is convenient to take the state at the mid point 
y = as a reference state and define the following dimensionless quantities: 



V* = V/vq, y* = yv' /v , /| = f e vl/n , 



(3.8) 



P* = P/Pa, u* = u/v , T* = T/T Q , g* = g/v' v 



(3.9) 



"V^o) p * = p /po, q* = q/po^o, 



(3.10) 



where, as in Eqs. I|2. 37(1 - 1)2. 42|l . the subscript denotes quantities at y — 0. In particular, vo = (2To/m) 1 ^ 2 is the 
thermal velocity Vth at y — 0. The reduced quantity y* measures distance in units of a nominal mean free path, while 
g* measures the strength of the gravity field on a particle moving with the thermal velocity along a distance of the 
order of the mean free path. The choice of l/v' (which depends on a) instead of 1/uq (which is independent of a) as 
the time unit is suggested by a larger simplicity in the calculations stemming from the kinetic model. In any case, in 
Section llVI we will summarize the results in real units, so the final expressions are independent of the specific choice 
of reduced quantities. 

In the above units, the kinetic equation l(3.7|l becomes 



fl + $) 



Vy'dy. log/,* 



T* 



,dul 



v dy 
T* d 
~2"dV' 



= \9 



V Q y * 



d 

W; 



dV 



(3.11) 
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where 



d\ogp* (V 2 5\ dlogT* 
dy* + yY ~ 2 J dy* 



d r log/; = + ( — - - ) -^r-. (3.12) 



On the right-hand side of Eq. (|3.11() we have taken into account that Q = Qv', where [cf. Eq. i|2.14[l ] 

, _ ^(1-a 2 ) 

^° — ar \ i 5~7i 2\ (3.13) 
0(a) + y^(l - a 2 ) 

is a pure number that only depends on the coefficient of restitution. It gives the cooling rate at any given point in 
units of the modified collision frequency v' at that same point. 

Our purpose is to solve Eq. I|3.11|l to second order in g* and get the associated hydrodynamic profiles. 

B. Perturbation expansion 

In this Subsection, all the quantities will be understood to be expressed in reduced units and the asterisks will be 
dropped. The expansion of <& in powers of g is 

$ = $ (1) .g + $ (2) .g 2 + C( ff 3 ), (3.14) 

where we have taken into account that the solution of Eq. 1)3. 5|) in the absence of gravity is / = fi with uniform n, u, 
and T. The expansions for the hydrodynamic fields have the forms 

p=\+pWg 2 + 0{g A ), (3.15) 
u z = u^g + 0(g 3 ), (3.16) 



T= 1+T (2) .g 2 + 0( 5 4 ). (3.17) 

Here we have taken into account that, because of the symmetry of the problem, p and T are even functions of g, 
while u z is an odd function. Also, without loss of generality, we have taken uq = 0, i.e., we are performing a Galilean 
change to a reference frame moving with the fluid at y = 0. Since v' = pT^ 1 / 2 , we have 

*/ = l+(V 2) - V)). 9 2 + 0( 5 4 ). (3.18) 

Nevertheless, only v 1 = 1 is needed in the evaluation of $W and $' 2 '. 

In order to solve Eq. 1)3.11(1 at each order, we will need to use the consistency conditions 

dVf e $ = 0, (3.19) 



dV V y fe$ - 0, (3.20) 



dVV z fe$ = 0, (3.21) 



dVV 2 fi<f> = 0. (3.22) 
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1. First order- 
To first order in g, Eq. I|3.11[) yields 



(1 - A) *« = -—^—V z ( 1 + V y ^ ) ee 0« (3.23) 



Wo* Z \ " dy 

where A is the operator 



2(1 - Co*) V 2 <9V / <9V 1-Co 



The function (j)^ has a known velocity dependence. Its space dependence occurs through ijW ; which remains unknown 
so far. In order to solve Eq. I|3.23|) . we will follow a heuristic procedure. First, we guess that the first-order velocity 
profile is parabolic: 

u«(l/)=4V. (3-25) 

Next, we note that the formal solution to Eq. (|3.23() is (f>( 1 ^ = J2kLo •A k (f>^ an d that the functional structure of 
A k (j)^ 1 ' remains the same for any k. Consequently, the solution to Eq. (|3.23ll must have such a structure, namely 

$ (1) (y,V) = V z (ao + ai V* + a 2 V y y). (3.26) 

Equations (|3.25|l and (|3.26|) have the same structure as the solution of the BGK equation in the elastic case^i^ 
Insertion of Eq. 13.26|) into Eq. (|3.23|) allows one to identify the coefficients ao, a\, a^. The result is 

«o = 4 So 2 -g , ai = — a 2 = -iu { ^. (3.27) 

4 - ko z + SO 

The consistency conditions H3.19|l . H3.20|l . and (|3.22|) are verified by symmetry. The coefficient u 2 is determined by 
the condition l|3.21[) with the result 

4 X) = I- (3.28) 
Once we know explicitly, we can get the non-zero components of the fluxes to first order. They are 

P$(y) = 2 J dVV v V z f $W = _ 2y (3.29) 
qV>{y) = J dVV'VJo^ = 2^r, (3.30) 

where f = n~ 3 / 2 e~ y2 is J t dXy = 0. 

2. Second order 

We proceed in a similar way as before. The equation for <j>( 2 ) is 

dvF>\ ( d 



S.0 



V 



dy \ 2 J dy 



(3.31) 



Wo 

Now, we guess the profiles 

P< 2 >(W)=*4V. (3-32) 
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r(%) = r 2 (V + riV. 



(3.33) 



The structure of A k (jr- 2 ' suggests the trial function 

^ {y, V) = b + b x V 2 + b 2 V y y + b 3 y 2 + b 4 V y 4 + b 5 V 3 y + b 6 V 2 y 2 
+b 7 V y y 3 + V 2 (co + ciV 2 + c 2 V y y + c 3 y 2 + Ci Vy 
+c 5 V 3 y + c 6 V 2 y 2 ) + V 2 (d Q + d.V 2 + d 2 V y y + d 3 y 2 
+d i V v i + d 5 V 3 y + d e V 2 y 2 + d 7 V y y 3 ) . 



(3.34) 



______ (2} (*2} f2^ 

Insertion into Eq. I|3.31|l allows one to get the coefficients bi, Ci, and di in terms of p 2 , T 2 , and T 4 . Condition 
(|3.2UI) is identically satisfied regardless of the values of p 2 2 \ and T^ 2 \ while Eq. (|3.21|) is verified by symmetry. 
On the other hand, conditions (|3.19|) and 1)3. 22[) yield 



(2) 
P2 



24 



-(2) _ 
1 2 — 



4 38 + 43C * + 17C * 2 
25 (l + C *)(2 + Co*) 



T (2) _ 
1 4 — 



15 



(2 + Co)- 



(3.35) 



The expressions of the coefficients 6^, Ci, and di as functions of a are given in Appendix^] From <I>( 2 ) we can calculate 
the second order contributions to the fluxes: 



P« -P (2) +2 JdVV 2 f ^ 2) 



*2 



24 102 + 87Co + 13Co 

25 (l + Co)(2 + Co) 2 ' 



(3.36) 



P${y) 



,(2) 



2 / dVlf/o$ (2) 



32 82 + 67Co +8Co 
25(l + G)(2 + Co*) 2 



56 



(3.37) 



q y 2 \y)= J dVV 2 V 2 f^ = \y 3 . 



(3.38) 



IV. SUMMARY AND DISCUSSION 



A. Hydrodynamic profiles 



Let us summarize here the main results obtained from the kinetic model through second order in the gravity field. 
The hydrodynamic profiles are given by Eqs. ft7r5 H 5T7 ) . (|3~23j) . (|3~2~1) . (|3~32l) . (|3~33l) . and iJOSJ- Expressed in real 
units, they are 



p(y) = Po 



Gfrng 
5 \T 



0(g% 



(4.1) 



u z (y)=u + !^y 2 + O(g 3 ), 
2t? 



(4.2) 



2 2 

1- ^ 
1277oK r 



1 38 + 43C * + 17C<f fmg 
25 (l + Co*)(2 + Co) \T 



0(g 4 



(4.3) 



In Eqs. (|4.2I) and (|4.3|l . 770 and k = Kg are the NS transport coefficients (evaluated at the mid point y = 0) of the 
granular gas heated by the stochastic force. In the kinetic model, those transport coefficients are given by Eq. (|3.2|l . 
Note that the elimination of the collision frequencies v or v 1 = j3v + Ci in favor of the transport coefficients r\ and 
k allows us to circumvent the limitation inherent to BGK-like models of not giving the correct Prandtl number. In 
that way, Eqs. (|4.1|) - (|4.3|) can be expected to be close to the Boltzmann results, as happens in the elastic casei— i- 



i n 




Therefore, in what follows we will use for rj and k the Boltzmann expressions (f2.33|) . In Eq. 1)4.3(1 the parameter Q is 
given by Eq. (|3.13l) . where (3(a) can be freely chosen. Here we will take the choice (|3.3fl . which makes the NS shear 
viscosity of the kinetic model agree with that of the Boltzmann equation. 

Comparison of Eqs. H4.lfl - H4.3fl witn tne NS predictions, Eqs. (|2.34fi . I|2.38fl . and (|2.39l) shows that the latter provide 
an incomplete description to second order in g. According to the kinetic theory description, the pressure increases 
parabolically from the mid layer rather than being uniform and the temperature has an extra positive quadratic term 
that is responsible for the fact that the temperature has a local minimum at y = rather than a maximum. This 
minimum is surrounded by two symmetric maxima located at a distance (j/ m ax) from y = of the order of a few mean 
free paths. Analogously, the NS equation of state 12.4Ufl is corrected by an extra term, 

rp rp Vo , , 2 , 1 38 + 43C * + 17Cp* 2 , , xn ,4, lA A , 

T = Tq -^-^ ] + 30 (1 + Co*)(2 + Co) {p - p ° )+0( * g } " (44) 

In order to analyze in detail the temperature profile (|4.3fl . let us measure the coordinate y in units of the mean 
free path Ao = (irV^nocr 2 )^ 1 — ^/by/ir^o/vo, where vq — -\/2To/m is the thermal velocity and vq is the collision 
frequency (|2.15fl . both at y = 0. Thus, Eq. I|4.4fi becomes 

where 

„ / x 4 , w ,„ „„ , , 4 2719 - 2741a + 706a 2 ,„ „, 

M^^^i^afiS-a^-SSa), ^25 (7-4a)(23-lla) ' (46) 

The coefficient ^(a) monotonically decreases with increasing inelasticity, while A^ia) has a maximum at a ~ 0.46, 
essentially due to the non-monotonic behavior of the thermal conductivity. 
The location y max of the two symmetric maxima is 



Note that, in the regime g\o /v 2 -C 1, y m ax is independent of the precise value of g. The relative value of the maximum 
temperature is 

T _ Al{a) (9±A 2 + 0{gA) (48) 



T 4 A 4 (a) V wg 



1 Q 




a 

FIG. 3: Plot of |y m ax|/Ao and (T max — Tq)/Tq versus a. In the latter case we have taken gXo/vo = 0.05. 



Of course, if we formally make Aa(a) = in Eq. 14.5(1 . the NS temperature profile 1(2.39(1 is recovered. 

As an illustration of the corrections over the NS description provided by the kinetic model, let us consider a 
value gXo/vQ = 0.05. In the case of terrestrial gravity, the above value corresponds, for instance, to Ao ~ 5 mm 
and vq ~ 1 m/s. Although terms of order higher than g 2 in 1(4.5(1 might not be negligible for this particular value of 
g\o/ v oi the qualitative features are expected to remain correct. Figure [21 shows the temperature profiles for a granular 
gas with a = 0.5 and a = 0.8, as well as for a gas of elastic particles (a = 1), as predicted by the NS and kinetic 
theory descriptions. We observe that strong deviations from the NS profiles are apparent, both for elastic and inelastic 
systems. Focusing now on the profiles predicted by the kinetic model, we se that, as the inelasticity increases, the 
locations of the two maxima shift towards the center of the slab and the value of the maximum temperature decreases. 
This behavior, however, is reversed if a < 0.4. The a-dependence of |y m ax|/Ao and (T max — Tq)/Tq is displayed in Fig. 
[3] The non-monotonic behaviors of y ma x and T max are consequences of that of A4 (a) . 



B. Fluxes 



The profiles for the elements of the pressure tensor and the components of the heat flux through second order are 
jiven by Eqs. ((3.29(1 . 1(3.30(1 . and 1(3.36(1 - 1(3.38(1 . Expressed in real units, they are 

P y M = -pagy + 0(g 3 ), (4.9) 
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P 



vv 



Po 



12 102 + 87C * + 13C *>o^S 



25 (1 + C *)(2 + G) 2 



Po 



(4.10) 



Pzz(y) = Po 



2^2 



16 82 + 67C * + 8Cr Po%9 



25(l + C *)(2 + Co) 2 



Pi 



14 



mg\ 2 



(4.11) 



2 2 

/ n Po9 3 

= -5 — 2/ 

3% 



(4.12) 



g 2 = -mK 5 + C(g 3 ). 
5 



(4.13) 



The ra-element of the pressure tensor is P xx — 3p— P yy — P zz . Note that P yy is uniform, in agreement with the exact 



balance equation (|2.22|) . Likewise, it is easy to check that Eqs. Q4.2|) . I|4.9|l . and 14.12|l are consistent with the energy 
balance equation (|2.24|l . Moreover, since the density profile is known through second order [cf. Eqs. I|4.1|l and l|4.3|) ]. 
Eq. i|2.23[l can be used to get P yz through third order. 

The shear stress P yz agrees to second order in g with Newton's viscosity law (12.26(1 . However, the component q y 
of the heat flux parallel to the thermal gradient does not obey Fourier's law l(2.27|l (note that /x « in the heated 
state). In fact, from Eqs. 1)4. 3[) and (|4.12|l one can write an 



Qy 



d_ 



T 



-V 2 T 



0(g% 



(4.14) 



which shows that one needs to incorporate super-Burnett contributions to account for the relationship between the 
heat flux and the thermal gradients. The extra term on the right-hand side of Eq. (|4. is responsible for the counter- 
intuitive fact of q y having the same sign as dT /dy in the region < \y\ < |y m ax|, i-e., the temperature increases as one 
moves away from the mid layer y — and yet the heat flows outward from the colder to the hotter layers. A steady 
state is still possible because the energy deficit is compensated for by the viscous heating. An additional departure 
from Fourier's law is related to the existence of a component q z of the heat flux normal to the thermal gradient, an 
effect that is already of first order in g and is related to a Burnett contribution associated with V 2 u z & 

Equations i|4.1[l . I|4.10|l . and (|4.11|) show that normal stress differences appear to order g 2 . It is easy to check that 
P yy < P xx < p < P zz , i.e., normal stresses are maximal along the flow direction and minimal along the direction 
normal to the plates. In order to characterize the normal stress differences, let us define the viscometric quantities 



p(y) 



p(y) 



(4.15) 



Their expressions are 



Mi/) 



1507T 



827- 733a + 158a 2 



(1 + a) 2 (23 - lla) 2 (3 - a)(7 - 4a) 



+8 



0(g% 



(4.16) 



Mi/) 



6?r 



38467 - 34763a + 7708a 2 
(1 + a) 2 (23 - lla) 2 (3 - a)(7 - 4a) 



56 



.9A0 



0(g 4 ). 



(4.17) 



Figure^shows the profiles of Ai(j/) and A 2 (y) for a = 0.5, a = 0.8, and a = 1 in the case gXo/v^ = 0.05. We observe 
that the normal stress differences increase with the separation from the mid layer y = 0. Moreover, those differences 
are more important for elastic gases (a = 1) than for inelastic gases (a — 0.8 and a — 0.5). However, as in the case 
of the quantities plotted in Fig. El the a-dependence of Ai and A 2 is not monotonic. This is illustrated in Fig. [3] for 
the point y = 0. We observe that the minimum values of Ai(0) and A 2 (0) occur at a w 0.5. 
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V. CONCLUDING REMARKS 



In this paper we have carried out a kinetic theory study of the steady planar Poiseuille flow undergone by a dilute 
granular gas under the action of the acceleration of gravity. In order to compensate locally for the energy loss due 
to the inelasticity of collisions, an external energy input in the form of a white noise driving has been assumed. This 
type of driving mechanism has been introduced in the literature to mimic the heating effects due vibrating boundaries 
without the complications associated with boundary effects. This is especially convenient in our approach, since we 
have been mainly interested in the bulk properties of the gas, namely in a slab centered in the middle layer having a 
width of the order of several mean free paths, away from the walls. 

Since granular gases are made of mesoscopic particles, terrestrial gravity (g — 9.8 m/s 2 ) plays in general a relevant 
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role, in contrast to the case of molecular gases. The dimensionless parameter characterizing the influence of gravity 
during the free flight of a particle between two successive collisions is g\/v^ h , where A is the mean free path and v t h 
is a typical (thermal) velocity. Under many conditions of practical interest, the parameter g\/v 2 h can have a non- 
negligible effect and yet be sufficiently small as to justify a perturbative treatment. For instance, if A ~ 1 mm-1 cm and 
^th si 1 m /s, which are typical values in experiments on metallic or glass spheres, one can have g\/v 2 h ~ 10 -3 -10 -1 . 
Therefore, in our study we have performed a perturbation expansion of the velocity distribution function in powers of 
the gravity strength through second order. The reference state (i.e., the state at zero gravity) is the steady uniform 
state heated by a white noise thermostat. Since the Boltzmann equation for inelastic spheres is quite complicated 
to deal with, we have employed a kinetic model equation inspired in the BGK model. This has allowed us to 
obtain explicitly the velocity distribution function through second order in terms of the velocity vector, the spatial 
coordinate, and the coefficient of restitution. By velocity integration one can obtain any desired moment, but here 
we have focused on the hydrodynamic fields (pressure, flow velocity, and granular temperature) and their associated 
fluxes (stress tensor and heat flux vector). 

The results show that the non-Newtonian features previously studied in the case of elastic particles^S*2iiiiiii4ii& 
persist when inelasticity is present. In particular, the temperature profile T(y) exhibits a bimodal shape: it has a local 
minimum To at the central layer and reaches two symmetric maxima T max at a distance |?/max| of about three mean 
free paths. The relative height of the two maxima, (T max — Tq)/Tq is about 10 times the square of the dimensionless 
parameter gX/v 2 h . On the other hand, the heat flows outward from the central layer, so it goes from the colder to 
the hotter layers within the region \y\ < \y max \. Other non-Newtonian effects include normal stress differences and 
the existence of a component of the heat flux parallel to the flow and hence normal to the thermal gradient. 

The fact that the nonlinear transport properties of the granular Poiseuille flow are qualitatively similar to those of 
the elastic case does not come as a surprise, especially since the characteristic collisional cooling of the granular gas 
is balanced by an external driving. In that context, our aim in the present work has been two-fold. On the one hand, 
the example of gravity-driven Poiseuille flow allows one to emphasize once more that granular gases constitute an 
excellent playground to reveal interesting (and even counter-intuitive) non-Newtonian phenomena on scales accessible 
to laboratory conditions. More importantly, we wanted to assess the influence of inelasticity on the departure of 
the Poiseuille profiles from the Navicr Stokes predictions. This influence is not easy to foretell a priori by means of 
intuitive or hand-waving arguments. According to the results reported in this paper, for small or moderate inelasticity 
(say a > 0.5) there is a slight decrease in the quantitative deviations from the Navier-Stokes profiles as inelasticity 
grows: the two temperature maxima becomes lower and closer, while the normal stress differences become smaller. 
The opposite behavior takes place for high inelasticity (a < 0.5), although that range is less interesting from an 
experimental point of view. 
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APPENDIX A: EXPRESSIONS FOR THE COEFFICIENTS bi, a, AND d t 

In this Appendix we list the explicit expressions of the coefficients in the expression for the velocity distribution 
function to order g 2 , Eq. I|3.34fl . They are 

h 4 (2 + 5Co*)(5 + 7Cq*+6C * 2 ) 

5(l + Co*) 2 (2 + Co)(l + 2Co)(2 + 3C *r 1 ' 

_ 8 160 + 622C * - 1051Co 2 - 2829Cq 3 - 1696<Co 4 - 276 C * 5 f A2l 

1_ 25 (l + Co*) 2 (2 + Co) 2 (2 + 7Co*+6C * 2 ) ' 

48 1 -h 3£* 

fc = -T(i + Co*)(2 + C5)(2 + 3C5)' ' 3 = °' (A3) 

b 32 5 + 29Co* + 12C * 2 ) 

4 5 (l + Co*)(2 + Co*)(l + 2Co*)(2 + 3C *)' 1 ' 
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32 (5 + 2G)(l + 3G) _ 16 Co* , = _8 

5 5 (l + a)(2 + Co*)(2 + 3C *)' 6 5 1 + Co*' 7 3' 



_ 4 + 12Cp* + 113Cq* 2 + 176Cq* 3 + 79Co* 4 + 6Q 5 . 
(l + Co*) 2 (2 + Co*) 2 (l + 2a)(2 + 3C *) ' 1 ' 



(3 + 2C *)(2-llC *~12C * 2 ) 
(l + a) 2 (2 + Co) 2 (l + 2Co*)(2 + 3C *)' 1 ' 



8-22C *-23Co* 2 -3C * 3 _ g 

C2 ib (l + C *)(2 + C *) 2 (2 + 3C *)' C3 8 1 + C *' 1 ' 



3 + 2C* 

C4 = 64 (1 + Co)(2 + Co)(l + 2Co)(2 + 3C5) ' (A9) 



3 + 2C* 16 

C5 = _64 (i + Co)(2 + Co)(2 + 3a)' C6 = iTa' (A10) 



, _ 8 76 - 48C * - 137C * 2 + 7C * 3 + 42( * 4 

25^°(l + a) 2 (2 + Co) 2 (l + 2Co)(2 + 3C *)' 1 ' 



_ 16 76 - 48g ~ 137C * 2 + 7( * 3 + 42C * 4 

1 25(i+c *) 2 (2 + a) 2 (i + 2a)(2+3a)' 1 1 



16 76 + 40C * + 23C * 2 + 21Co* 3 rA1 - 
2 25 (l + Co)(2 + Co*) 2 (2 + 3C *)' 3 5 1 + Co*' { ' 



64 1 



5 (l + Co*)(l + 2Co*)(2 + 3C *) 



* Electronic address: mtij@fsmek. ac.ma| 

T F 1 ,lpr*tT-nnir aHHrpw anrlrpti(Q)iinPY pq TT 



(A14) 



a - 64 1 , 16 1 . 16 

rf5 -T(i + C *)(2 + 3C *)' 6 = "TTTCo*' 7 "i^- (A15) 
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